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Abstract. We solve analytically the Kadanoff-Baym equations for a noninteracting junction 
connected to an arbitrary number of noninteracting wide-band terminals. The initial equilibrium 
state is properly described by the addition of an imaginary track to the time contour. From the 
solution we obtain the time-dependent electron densities and currents within the junction. The 
final results are analytic expressions as a function of time, and therefore no time propagation is 
needed — either in transient or in steady-state regimes. We further present and discuss some 
applications of the obtained formulae. 



1. Introduction 

The Landauer-Biittiker formula P, 0| provides an intuitive physical picture of the steady-state 
current flowing in a multi-terminal junction and it is simple to implement. First one calculates 
the steady-state current I a /3 in terminal j3 carried by the scattering states originating from 
terminal a ^ f3 and populated according to the electrochemical potential [i a . Then one sums 
the difference I a p — Ip a between the currents flowing in and out terminal (3 over all terminals 
a (3. This gives the steady-state current Ip in terminal (3. 

The first microscopic derivation (based on the time-dependent Schrodinger equation) of the 
Landauer-Biittiker formula was given by Caroli and co-workers 0, 0|. They considered the 
terminals initially uncontacted and in equilibrium at different chemical potentials. Then they 
switched on the contacts and derived the Landauer-Biittiker formula as the long-time limit of 
the expectation value at time t of the current operator. We will refer to this procedure as the 
partitioned approach. 

An alternative approach, more akin to the the way the experiments are carried out, was 
proposed by Cini about a decade later He considered the system initially contacted and 
in equilibrium at a unique chemical potential and then drove the system out of equilibrium by 
applying a bias voltage between the terminals. We will refer to this procedure as the partion-free 
approach. In both approaches one recovers the Landauer-Biittiker formula due to the loss of 
memory of the initial preparation 0| . 



The microscopic derivation of the Landauer-Biittiker formula requires the evaluation of the 
expectation value Ip(t) = Jg|\l/(i)) where is the many-body state of the system at 

time t and Ig is the current operator. Since the electrons are noninteracting this expectation 
value can be rewritten as the sum over all occupied one-particle states \ijj^{t)) = e~ lht \ipk) of 
Here h is the Hamiltonian of the contacted and biased system whereas \tpk) are 
the eigenstates of the Hamiltonian ho which describes either the non-biased uncontacted system 
(in the partitioned approach) or the non-biased contacted system (in the partition- free approach). 
For the evaluation of (V ; fc(^)|^|V ; fc(^)) one could naively insert a complete set of eigenstates \<f) q ) 
of h and evaluate the overlaps {ijjk\(j) q ) . This procedure is, however, numerically lengthy and 
unstable due to the singular 5-like contribution to the overlaps. The calculation of Ip(t) is 
most easily carried out using nonequilibrium Green's functions 0, @|- This mathematical tool 
when applied to quantum transport in multi-terminal junctions provides a natural framework to 
calculate the current at all times and not only at the steady state. 

In fact, there have been several attempts to generalize the Landauer-Biittiker formula to the 
time domain. Here we mention the work of Pastawski who derived a formula for Ip (t) using the 
partitioned approach in the linear response and adiabatic regime 0| . An important step forward 
in the calculation of Ip (t) was done by Jauho et al. |8j. These authors used the partitioned 
approach to write Ip (t) as a double integral (over time and energy) of the trace over the junction 
degrees of freedom of a calculable combination of Green's functions in the same region. In the 
special case of terminals with a wide band and of junctions with one single level it is possible to 
perform the time-integral and obtain a time-dependent version of the Landauer-Biittiker formula. 
This formula was then derived in Ref. 0| using the partition-free approach, thus confirming the 
loss of memory of the initial preparation. 

The derivation of a time-dependent Landauer-Biittiker formula for arbitrary junctions would 
be extremely useful to interpret the oscillations and damping times typically observed in the 
transient current after the sudden switch on of a bias. A progress in this direction was done in 
Ref. [l~C)| where the authors derived a time-dependent Landauer-Biittiker formula for the spin 
current of a single-level junction. 

In this work we generalize the results of Ref. [10] to junctions of any shape and dimensions 
using the wide-band limit approximation (WBLA) for noninteracting electrons (Sees. [2] and 
[3]). Furthermore we also derive a general formula for the time-dependent one-particle density 
matrix which can be used to calculate the local density and current density. We will work in 
the partition-free approach which is conceptually easier since it does not involve the subtle issue 
of different chemical potentials in equilibrium. The final formulae for the current and the one- 
particle density matrix have the merit of elucidating the relative importance of the electronic 
transitions at a certain time. As an illustration we will use these formulae to calculate the 
transient response of a ring-shaped junction (Sec. HJ. 



2. Assumptions and set-up 

We investigate the following quantum transport setup: An arbitrary number of metallic leads 
(a) acting as charge-carrier reservoirs are connected to a lattice network acting as a molecular 
device (C). We assume that the electron transport is ballistic and therefore neglect the electron- 
electron interactions. We will also assume that the energy eigenvalues of the Hamiltonian of the 
molecular device are well inside the continuous energy spectrum of the leads and use the WBLA. 
The described set-up is characterized by the following Hamiltonian: 

H = ^2 e ka™ka,a + ^ ^rnn^rn,a^n,a + ^ \^mkJ^ra^ka,(j + ^kam^ka^m^ ■ (1) 
ka,a mn,cr mka,cr 



The first term accounts for the leads with ka indexing the k:th. basis function of the a = 1,2,3,. 



lead. The single-particle spectrum of the leads is e ka and the number operator in the leads is 
expressed in terms of the creation and annihilation operators as n 



fcay 



d\ d, , with a the 

spin index. The second term is for the molecular device, or central region, (indices m and n) 
with creation and annihilation operators dXi,<T and d m ,a and hoppings T mn between sites m and 
n. The last term is for the coupling between the central region and the leads with hoppings 

At times t < to the system is in thermal equilibrium at inverse temperature /3 and chemical 
potential //, the density matrix having the form p = ^q~^\ h ~i iN ) where Z is the grand-canonical 
partition function. At t = to the lead energy levels are suddenly shifted by some constant value, 
£fca —> tka + V a , to model the sudden switch-on of an external bias voltage in the a:th lead. 
This means that the system is driven out of equilibrium and charge carriers start to flow through 
the central region. To calculate the time-dependent current we use the equations of motion for 
the one-particle Green's function on the Keldysh contour 7k- This quantity is defined as the 
ensemble average of the contour-ordered product of particle creation and annihilation operators 
in the Heisenberg picture 

G rs (z,z') = -i(r 7K K iH (^)4 jH (z / )]) (2) 

where the indices r, s can be either indices in the leads or in the central region and the variables 
z, z' run on the contoui0. The matrix G with matrix elements G rs satisfies the equations of 
motion 



dz 



G(z,z>) 



h(z) 



dz' 



G(z, z) = 5(z, z)l , 



h{z>) 



S(z,z')l , 



(3) 



(4) 



with Kubo-Martin-Schwinger (KMS) boundary conditions. Here h(z) is the single-particle 
Hamiltonian. In the basis ka and m the matrix h has the following block structure 



(5) 



where {h aa i)kk' = baa'&kk'tka corresponds to the leads, (h a c)km = T^kam is the coupling part, and 
(hcc)mn = T mn accounts for the central region. As the system is initially in thermal equilibrium 
we have that for z on the vertical track of the contour €k a (z) = eta — Hi Tfc Qm (z) = Tk am 
and T mn (z) = T mn — p5 mn . On the other hand for z on the horizontal branches we have 
£ka(z) = eta + Va, T kam (z) = T kam and T mn (z) = T mn . Due to the coupling between the central 
region and the leads the matrix G has nonvanishing entries everywhere 
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(6) 



In the next Section we solve the equations of motion ([3|) and (jl]) for the Green's function Gcc 
projected in the central region. 

The contour has a forward and a backward branch on the real-time axis, [to, oo[, and also a vertical branch on 
the imaginary axis, [to, to — i/3] with inverse temperature /3, see e.g. [111 ). 



3. Derivation of the time-dependent density and current 

3.1. Projecting the equation of motion 

We project the equation of motion ([3|) onto regions CC and aC. The equation for G a c can 
be integrated using the Green's function g aa (z,z') of the isolated a:th reservoir. This Green's 
function solves the equation of motion [i^ — h aa \ g a a(z,z') = 5(z,z') as well as the adjoint 
equation with KMS boundary conditions. Introducing the embedding self-energy (with indices 
in region C) 

£em(z,z) = ^E a (z,z') ; E a (z, z') = hca9aa(z, z')h a C (7) 



we obtain the equation of motion for the Green's function of the central region (the subscripts 
CC are omitted from now on) 



dz 



h 



G(z,z') = 6(z,z')+ / dzS em (z,z)G(z,z') 

J 7k 



The adjoint equation of motion can be derived similarly and read ll| 



G(z,z') 



dz' 



h 



S(z,z') + 



/ dzG(z, z)S e 



\z,z 



(8) 



(9) 



The embedded equations of motion for G have the same structure as the Kadanoff-Baym 
equations (KBE), the difference being that the many-body self-energy is replaced by the 
embedding self-energy. In the case of interacting electrons with an interaction only in the central 
region Eqs. ([8]) and Q are modified by the addition of the many-body self-energy S to the 
embedding self-energy S em , i.e., I7 em — > E em -\- S. Since S = S[G] is a functional of the Green's 
function in region C the embedded equations of motion in the interacting case constitute a 
closed set of integro-differential, nonlinear equations for G The simplification brought by 

the absence of interactions is that the KBE ([8]) and Q are linear in G since the embedding 
self-energy is completely specified by the parameters of the Hamiltonian. 

The density and current density can be extracted from the lesser component of the Green's 
function at equal time. We denote by z = i_ the contour point on the forward branch, z = t + 
the contour point on the backward branch and z = to — ir the contour point on the vertical track. 
The Keldysh components lesser (<), greater (>), retarded (R), advanced (A), left (J), right (]) 
and Matsubara (M) of a function k(z, z') on the contour are defined according to [12] 
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(10) 

(11) 

(12) 
(13) 
(14) 
(15) 
(16) 



To generate an equation for G < we subtract Eq. ([9]) from Eq. ([8]) and set z = z' = t',. Taking 
into account that h(z) = h is independent of z for z on the horizontal branches we get at equal 
time 



i^G<(t,t)- [h,G<(t,t) 



[S em ■ G 



G<-St] (t,t) 



+ 



(t,t) , 



(17) 



where we defined [f-g](t,t') = /~ dif(t,t)g(i,t') and [f*g](t,t) = -itf drf(t,T)g(T,t'). 
Equation (11711 can also be written as 

i—G<(t,t)- [h,G < {t,t)] = - [g r • S< m + G< ■ + GU 4m] (t,t)+h.c. (18) 

where we used the properties of G and S era under complex conjugation 

Let us comment Eq. f)18|) briefly. Setting the right-hand side to zero we see that Eq. (I18p 
reduces to the Liouville equation for the one-particle density matrix p = — iG of the isolated 
central region. Thus the embedding self-energy accounts for the openness of region C. The first 
term inside the square brackets is a convolution between the propagator in region C, G R , and 
Sf m . Since Sf m is proportional to the probability of finding an electron in the leads this term can 
be interpreted as a source term, i.e., a term that describes the injection of electrons into region 
C. The second term has the opposite structure: a propagator in the leads, S em , is convoluted 
with G < which is propotional to the probability of finding an electron in region C. Thus this 
term can be interpreted as a drain term and is responsible for damping and equilibration effects. 
The last term inside the square brackets accounts for the initial preparation of the system. In the 
partioned approach this term would be zero since the hopping integrals T^ am = in equilibrium. 
However, in the partition-free approach this term is nonzero and accounts for the initial coupling 
of the central region to the leads. 

More generally convolutions along the vertical track carry information on the initial 
preparation of the system. For instance for a system of interacting electrons we can either start 
with a noninteracting system and then switch on the interaction in real time or we can start 
with a system already interacting. In the latter case the many-body self-energy is nonvanishing 
on the vertical track and the convolution G * E accounts for the effects of initial correlations. 

3.2. Self- energy and Green's function calculations 

The solution of Eq. (fT8|) requires first to calculate the Matsubara component G M , and then 
from G M the right and left component G"l and G^ . The Matsubara component G M can be 
determined from the retarded/advanced components by analytic continuation, see below. Since 
the equations for G M , G^ and G^ contain the embedding self-energy the preliminary step is to 
obtain an expression for S em . 

Having a time-independent Hamiltonian (on the horizontal branches of the contour) the 
retarded/advanced components of the self-energy depend only on the time difference 

k 

where g^ a is the diagonal element of the Green's function g aa of the isolated a:th lead, see Eq. 
([7]). The retarded component of the self-energy is found by conjug ating Z*(t,t') = [S^(t',t)]\ 
According to the WBLA the Fourier transform of £^ i s frequency independent 

^a.mnH = ^ ] -^mfca 9ka (^0 ^kan = ^ ] ^mfca 77 '■ ^fcan = 7^^ a,mn (^0) 

k k 

which implies that i s a ^ so independent of the external bias voltage V a . The time-dependent 
self-energy of Eq. (|19p is therefore 

<mn(M') = / ^e-^'^Ju;) = ^r a , mn S(t - . (21) 



Within the WBLA we can easily calculate the two other self-energy components in Eq. (fl~8|) (see 

do; e^ 4 ^)* 



Appendix A) 



^ a,mn (^~: ^) 



r 

- 1 a,rr 
i-^ry.? 



2ir 0J q — UJ + {J, 



2^ 



(22) 
(23) 



where the sum over q is a sum over the Matsubara frequencies oj q = ^ffl* , and the function / 

is the Fermi function, f(u) = l/[e^ w + 1]. 

Having the explicit form of the self-energy components we can derive expressions for 
the Green's function. For the following calculations it is convenient to define the effective 
Hamiltonian h e g = h — =4> = h+ ^T, where r = ^ a r a . This effective Hamiltonian 
is therefore non-hermitean. The two Green's function components in the square brackets of 
Eq. ([18]) read (see |Appendix B[ ) 



G\t t r) 
G K (t,t') 



-ih cB t 







G M (0,r)-/ dt'e ih ^' dfS] m (t',f)G M (f,T 



-W(t - t')e 



-ih eS (t-t') 



(24) 
(25) 



with G M the Matsubara Green's function. G M can be obtained from G R and G A by analytic 
continuation since G M (oj q ) = G R (uj q + fi) if Imfc^] > and G M (oo q ) = G A (u q + /i) if Im[u;q] < 0. 
see |Appendix B" 



Now we have all ingredients to calculate the convolutions in Eq. (|18p . We report here the 
final results and refer to Appendix G for details. The three terms read 



[G R -S< m ](t,t) = i£/i±i/( w - M ) 



duj 



1 -e 



i(u}+V a -h cB ) 



G R (u; + V a )r a , 



[G<-stm](t,t) = ~G<(t,t)r 



dw 



GUEL (t,t) = i/^/^-riEe^-^H^ 



(26) 

(27) 
(28) 



3.3. Solving Eq. (EJ for G<(t,t) 

We insert Eqs. (|2"fi ft . (|2T )l and pH ]) into Eq. <^8§ and get 



ij t G<(t,t)-[h,G<(t,t)} 



^i(u}+V a -h eS )t 



G R (oj + v a )r a + ^-rG < (t,t) 



+ - M) E e^- h ^G^)r a | + h.c. 



(29) 



This is a nonhomogeneous, linear, first-order differential equation for G K (t, t) and, therefore, can 
be solved explicitly. The solution is worked out in Appendix D and reads 



iG<(t,t) = [^f(u-n)Y,{ A a {u + V a ) 

J a 

' e i(^+V a -h eS )t G R^- )Aa ^ + + h c 



+ V a 



+ V^- iKst G K (u)A a (u + V a )G A {uj)e<^ } , 
where we introduced the partial spectral function as 

Aa(u) = G R (u)r a G A (u) . 



(30) 
(31) 



The full nonequilibrium spectral function is A{oS) = ^ a A a {u). 

Given the original complexity of the problem the final result is surprisingly compact. Equation 
(|30|) is an explicit closed formula for the equal-time G < or, equivalently, for the reduced one- 
particle density matrix. All the quantities inside the integral can be calculated separately, and 
no time-propagation nor self-consistency algorithms are needed. Also, we may extract several 
physical properties: 

(i) With no external bias, V a = 0, only the first row contributes. This term correctly gives 
the equilibrium value of the equal-time G < since at zero bias ^2 a A a (uj) is the equilibrium 
spectral function. 

(ii) Both the second and the third row vanish exponentially in the long-time limit, and the 
equal-time G < approaches a unique steady-state value. 

(iii) The transient dynamics is given by the second and the third row. By inserting a complete 
set of eigenstates of the effective Hamiltonian h e s we notice that: 

(a) The second row gives rise to oscillations with frequency coj = \fi + V a — e^ ff | where 

is the real part of the j:th complex eigenvalue of h e Q. These oscillations correspond to 
transitions between the biased Fermi level of the leads and the resonant levels of the 
central molecule. 

(b) The third term accounts for intramolecular transitions and leads to oscillations with 



frequency uj 



. These oscillations are visible only if the effective Hamiltonian 



h e s does not commute with r a . In the case that [h e ff, r a ] = the time dependence of 
the third term is of the form e _lfeeR * +lfe efi* = e - -^'. 

3.4- Current calculation 

The time-dependent current through the interface between the central region and the a:th 
reservoir is calculated from the following equation fill ] 



I a (t) =4( Z Re{Tr \s< ■ G A + ■ G< + S\ * gH (t,t)\ . 



(32) 



For the terms inside Eq. (|32|) we proceed in the same manner as we did previously to obtain the 
results in Eqs. (EgJ, |E2J and PHI): 



[£<-G A ](t,t) = il^-f(u>-»)r a G A (oo + V a 



1 



[E*'G<] (t,t) 



--r a G<(t,t) 



El*Gt]{t,t) = i / ^f(u - /x)r a G A (o;)e- i ^-^ 



(33) 
(34) 
(35) 



Figure 1: Six-site ring coupled symmetrically to one-dimensional TB leads. 




Figure 2: Schematic of the transport setup through the energy levels of the molecule. 



Inserting these results into Eq. (|32p and taking into account the explicit expression for G < (t, t) 
in Eq. flU we get 



r a G K (u + v p )r p G A {uj + Vp) - r a G K (oj + v a )r p G A (cj + v a ) 



+ V B 



r a j^ +v P- h ^ t G K {u) (-i6 a pG R (u + V p ) + A p (u + V p )) + h.c. 



+ vjr a e 



Ah 



■ st G K (uj)A f} (uj + Vf } )G A (u;)e ih ls t } 



(36) 



The physical interpretation of the terms in Eq. (|36|) is similar to the one after Eq. (|30j) . We have 
a steady-state part given by the first row, and a time-dependent part given by the second and 
the third rows. The time-dependent part vanishes exponentially in the long-time limit and the 
oscillations in the current have the same structure as in the reduced one-particle density matrix. 



4. Results 

Let us consider a six-site tight-binding ring connected to two tight-binding, semi-infinite, one- 
dimensional leads as shown schematically in Figs. [1] and [2j The parameters according to the 
notation in the figures are tc = —2.0 (hopping in the molecule), chemical potential \i = and 
zero temperature (/3 — > 00). We choose the hopping ti = ijj in the left/right lead to be much 
larger than any other energy scale. Then T, A (u ~ pt) = it^ c /t a + 0(l/t^) where t a c is the 
hopping between the molecule and the leads. For this situation the WBLA with r a = 2t^ c /t a 
is a very good approximation. We study the weak coupling case Tl^r = 0.1 and drive the system 
out of equilibrium by the sudden switch-on of a bias Vl = V = —Vr. We analyze the contribution 
of different terms in the charge current corresponding to different physical features as discussed 
below Eq. (|3U|) . In Eq. (|36|) the first row is 'steady state', the second row consists of '1st term, 
a' and '1st term, b' and the third row is '2nd term'. The second row is divided into two parts 
[~ G R (w)G R (a; + Vr) and ~ G^(co)Ap(uj + Vr)\ since they give rise to different features. 

In Figs. [3] and 0] we plot the current through the right interface and see that weakly biased 
leads, V = 0.5, give a negligible steady-state current. Transitions between the biased leads and 
the molecule are captured by the '1st term, a'. This is confirmed by the peak in the Fourier 
spectrum at uij = e^ ff ± V. Transitions between the molecular levels are accounted for by the 
'2nd term', as it can be seen in the Fourier transform with a peak at oj = 6. In addition to our 



previous observations: (1) the '1st term, b' also gives rise to intramolecular transitions and (2) 
there seems to be no intramolecular transitions at uj = 2 or u = 4. By expanding Eq. (|36p in the 
eigenbasis of the effective Hamiltonian and manipulating the terms further one can show that 
the the '1st term, b' contains a term of the form e -1 ^ _€fc '* which explains the first finding. 
The second finding suggests that there is some underlying selection rule for some of the energy 
levels and hence that some levels do not participate to the transport process. 
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Figure 3: Different terms of the time- 
dependent current in units of V through the 
right interface with symmetric coupling and 
weak bias. 



Figure 4: Absolute value of the Fourier 
transform of the terms of the current in 
units of T. 



If we increase the bias window to cover the first molecular levels, V = 2.5, then we see in 
Fig. [5] that the current has a non-zero steady-state value. Similar findings, as with weaker bias, 
for the possible transitions are seen in Fig. [U] We also see that there is a small bump at uj = 4 in 
'1st term, b' and '2nd term'. Given that the setup is completely identical to the previous case, 
this fact is due to a second (or higher) order response since the same symmetry arguments apply. 




Figure 5: Different terms of the time- 
dependent current in units of V through the 
right interface with symmetric coupling and 
strong bias. 

If we, however, distort the symmetry of the junction then also the intramolecular transitions 
with lower energies become visible. This is clearly seen in Figs. [7] and [8] where we connect the 
molecule asymmetrically to the leads (1st site to the left and 3rd site to the right, see Fig. [1]). We 
can also break the symmetry by deforming the molecule with, for instance, one hopping (between 
sites 1 and 6) being 2tc- This splits the degenerate levels in Fig. [2] and also the corresponding 
intramolecular transitions can be seen in Figs. l9land[T0l 



Figure 6: Absolute value of the Fourier 
transform of the terms of the current in 
units of r. 
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Figure 7: Different terms of the time- 
dependent current in units of T through the 
right interface with asymmetric coupling 
and weak bias. 



Figure 8: Absolute value of the Fourier 
transform of the terms of the current in 
units of r. 



Contact: LI <-> 1, 4 <-> Rl, deformed 

\ V L = -V R = \t c \li 



— 1st term, a 

— 1 st term, b 

— 2nd term 

— Full result 
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Figure 9: Different components of the time- 
dependent current in units of V through the 
right interface of a deformed molecule with 
symmetric coupling and weak bias. 
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Figure 10: Absolute value of the Fourier 
transform of the terms of the current in 
units of T. 



As the contributions from different terms sum up to the total current we can plot the full 
results for, e.g., the right current of the symmetrically coupled molecule against, e.g., the bias 
or the coupling strength. In Figs. fTTT [T2l [T3l and [H] we display the full results. The transient 
dynamics is visualized better but distinguishing between the different contributions is more 
complicated. In Fig. [13] and [T31 the axes are not scaled due to varying T. It is clear that by 
increasing the bias window more levels open up for transport and therefore the steady-state 
current grows. The oscillation frequencies corresponding to transitions between molecular levels 
remain unchanged while the oscillation frequencies corresponding to transitions between the 
molecule and the leads vary (peak shift). By increasing T, and hence by widening the resonances, 
electrons can flow even with intermediate bias voltages. Correspondingly, the steady-state value 
of the current increases, the relaxation time decreases whereas the oscillation frequencies remain 
invariant. 

5. Conclusions 

In conclusion we solved the Kadanoff-Baym equations for the Green's function of an open 
noninteracting system by properly taking into account the initial contacts between the system 
and the reservoirs. We used the analytic solution for the time-dependent density matrix to 
derive a time-dependent version of the Landauer-Biittiker formula. As an application we 




Figure 11: Time-dependent current in 
units of r through the right interface with 
symmetric coupling (T = 0.1) and varying 
bias. (Dotted lines refer to steady-state 
values.) 
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Figure 12: Absolute value of the Fourier 
transformed right current in units of T. 
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Figure 13: Time-dependent current 
through the right interface with symmet- 
ric coupling (varying strength) and bias 
Vl = —Vr = 2.5. (Dotted lines refer to 
steady-state values.) 



Figure 14: Absolute value of the Fourier 
transformed right current. 



considered a tight-binding benzene-shaped junction and calculated the time-dependent current 
flowing through it. The advantages of having an explicit solution are that the numerical effort is 
drastically reduced and that the transient behavior can easily be interpreted in terms of virtual 
transitions and decay rates. Our time-dependent Landauer-Biittiker formula holds promise for 
studying the transient behavior of large junctions like, e.g., wide nanoribbons or large-diameter 
nanotubes, as well as disordered junctions where a large number of simulations is required to 
perform the average over different configurations. 
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Appendix A. Self-energy calculations 

The Matsubara self-energy is an antiperiodic function (we are studying fermions) with period 
given by the inverse temperature j3. For the calculation of the Fourier coefficients we can use the 



relation: S^{oj ± ir]) = S^ A (lo + /i). Therefore the Matsubara self-energy is simply given by 

Z%mn(T,T') = ^E e " U ' (TV) CK) ( AJ ) 

where 



f -±r a Im[u q ] > 
^UK) = < (A.2) 
[ +^r a Im[w,] < 

and w g = (2q + l)?r/( — are the Matsubara frequencies. For the isolated Green's function of 
the biased a:th reservoir we have 

gUt,t ) = -i^-^e-^+^-'o) , (A.3) 

^ (r ' r/) = -ifl^ J ; + „ ■ (A - 4) 

Without loss of generality we take the time to a t which the bias is switched on to be zero. Then 
we can write 

gl a (t,r) = 0)^(0, r)=e-^ +v ^gf a (0,r) , (A.5) 

gl a (r,t) = -igf a (r,0)gt(0,t)=e i ^ +v ^gf a (T,0) . (A.6) 

By using Eqs. (1A.5|) and (|A.6P we can calculate the right and left embedding self-energies 



1 g-i(efea + Va)t 

— J P ~T ; W <? ~~ £ fc" T M 

Q k ^ 

1 ^ , /"did e ^ +y -)* 



9 



J P 7 ITT 0O q — UJ + jJL 



— j p ~r i w <? ~~ + m 

'do; e i(aJ+Kl) * 



r a , mn — ^ e"^ r / — , A.8 

— lp J 2tt cj q — uj + 



where we used F amn = 27r T m k a 5{uj — e^a — ^a)Tkan- It only remains to calculate the lesser 
component. We have 

= */(<*« " ^e-^^-^ (A.9) 

and therefore 

= iTa^f^fiu-rie-**™*-*). (A.10) 



Appendix B. Green's function calculations 

The Fourier coefficient of the Matsubara Green's function reads 



g m K) 



UJq — h+2 



u q -h- S em (uj q ) + (i 

lm[uj q ] > 
Im[u q ] < 



Vp^ ImK]>0 

km Im N <0 



(B.l) 



where we defined hgg = h — ^T. The right component of the Green's function can be derived 
from the equation of motion 



dt 



poo r/3 

G\t,r)= diSf m (t,t)G\t,r)-i dfS] m (t,f)G M (f,r). (B.2) 
Jo Jo 



The insertion of the retarded self-energy from Eq. (|20p leads to 



G\t, 



t = e 



-i/l e fft 



G M (0,r)- fdt'e^' dfE] m (t',f)G M (f,r) 



(B.3) 



where we noticed that G M (0,r) = C*1(0,t). 

Finally the retarded Green's function in Fourier space reads 



G (ui) 



u-h- sf m (uj) u-h + ^r 



and Fourier transforming back in the time domain we recover Eq. ([25 



(B.4) 



Appendix C. The three terms in Eq. ([18]) 

For the first term we use Eqs. (|25p and (jA.lOp to obtain 

/>oo 

[G R -S< m ](t,t) = / dtG R (t,t)S< m (i,t) 
Jo 



1 - e 



i(u)+Va— h e f{) 



G R (oj + V a )r a . (Gl) 



Also the second term is readily calculated by using Eq. ([21 



[G< ■ S^] (t,t) = J dtG<(t,t)St(t,t) = ^G<(t,t)r . 



(C.2) 



The third term involves somewhat more trickery because of the rather complicated form of the 
right Green's function. Inserting the expressions from Eqs. ([B.3P and (|A.8[> we get 



G^ ~k ^7l~ 



rP 

(t,t) = -i / drG\t,r)si m (T,t) 
Jo 



-ih eS t 



G M * uL 



(0, t)-i dt'e 
'0 



'Jh eS t' 



(f,*U. 
(C.3) 



By using dre^ 9 1 V) T = f35 qq ' for the Matsubara frequencies we may manipulate Eq. (|C,3P 
further. Inserting Eqs. ()A,7p . ()A,8j) and (jB.ip we obtain for the double convolution 



(t',t) 



i/ dr(-i) / dfS] m (t',f)G M (f,r)sUr,t) 



duj f du/ v - 1 v - 



0J q - UJ + /i 



G M K) 



e i(u/+V Q ,)t 

w 9 - w' + fl 



(C.4) 



The integration with respect to uj can be done by closing the contour in lower-half plane (LHP) 
because of the exponential convergence factor, whereas the integration with respect to uj' can be 
done by closing the contour in the upper-half plane (UHP). However, the corresponding poles 
are located on different half planes, and this makes the double integral to vanish for every uj q . 
Hence 



(t',t) = 



(C.5) 



and in Eq. (|C,3P we are left with 



(0,t) 



-i f drG M (0,r)sUr,t) 
Jo 



duj 1 



2vr -i/3 



G M (uj g )e^" 



UJq - OJ + \X 
q H a 



(w+V a )t 



(C.6) 



where on the last line a convergence factor e^ 9 was added to account for correct limiting 
behaviour when t — > 0. The sum over Matsubara frequencies can be performed using the 
Luttinger-Ward trick [13l | and yields 



— y 

-i/3 ^ 



G M (uj q )e^ 

i/3 ^ UJ q — UJ + /i 



L 2tt j v > uj' -uj + 



G m {uj'-\5) G M (uj' + i5) 



Mr, ,/ 



uj 1 — uj + /u — i<5 cj' — w + /i + i5 



(C.7) 



By inserting Eq. ()C.7p into Eq. ()C.6P we get 



(0,t) 



du; /" du/ 



2vr / 2vr 



G M (u/-i<5) G M (u/ + i<5) 



uj 1 — uj + jU — i<S u/ — cj + + i<5 



(C8) 

Now the integral over uj can be done by closing the contour in the UHP. The first term in square 
brackets integrates to zero because of the pole in the LHP. The pole of the second term occurs 
at uj = uj' + /i + iS, and therefore 



G M * Si 

<5 0+ 

UJ 1 = UJ — f_L 



(0,i) 



du' 



i / f(uj')G M (uj' + iS) e^'+» +i5+v "» 

i / ^f(uj - fj,) g m (uj - M ) £ rj^^ 



i(w+Va) 



(C.9) 



This can be inserted into Eq. (|C.3|) to obtain 



G\ * E[ m (t, t) = i / ^/( W - /*) £ e^-^G R ( W )r c 



(CIO) 



(D.2) 



Appendix D. Derivation of Eq. (|3"0|) 

We first state some useful identities for the retarded/advanced Green's function to be used later: 
G R/A (u; + V a ) = G R/A (co) - V a G R/A (co)G R/A (u + V a ) , (D.l) 
which can be checked directly by using Eq. (|B.4p and its adjoint. From Eq. (|D. 1[) it follows that 
V^G r (oj)G a (cj)G r (uj + V^G^ + V a ) = [G r (lo) - G r (oj + V a )] [G a (oj) - G a (uj + V a )] 

Taking into account Eq. (jD.lj) in Eq. (|29p we get 

i^G<(t, t) - h eS G<(t, t) + G<(i, t)^ 

-i f jLf(u -fx)J2 {e i ^- h ^V a G R (a;)G R (L. + V a ) + G R (^ + V a )) r a 
+ i f p.f(u> - (j.) ]T r a {g a (uj + y Q )G A (^)y Q e- i ^- /l e ft 0* + G A( W + . (D.3) 

It is convenient to rewrite the Green's function as G < (t,t) = e~ lhcBt G < (t,t)e lh es t . In this way 
the left-hand side of Eq. (|D.3j) becomes 



dt 



dt v ; 



(D.4) 



Then the right-hand side of Eq. (|D.3j) can be multiplied from left by e lh-eff * and from right by 



e A cff* to give 



iie< (M ) 



-i / p-f(w -tiY. e '""' [° R (" + v «) r « - r -° A ^ + v a)] 
-i / ^/(w -M)E y « [g r (u,)G r (u; + ^r^+^U 

-e- i(tJ+v ^-^ ff) *r a G A (w + y Q )G A (a;)j . (D.5) 
Now we are ready to integrate both sides over t to obtain 



G < (t,t) -G<(0,0+) = 
=G < (0,0+)=G M (0,0+) 



G R (a;)G R (w + V a )r a e i ^ +v "- h ^ t ' 
_ e -i(a;+K 1 -h eff )t' ra G A (w + y a )G A (w)]} . (D.6) 



The integration over t' for the second term in Eq (|D.6|) can easily be done. For the first term we 
need the following result: Given two arbitrary matrices A and B 



dt'e [At ' 



o 



1 -B - B- 1 



x - A x- At 



Utf = _ ie iAt> _J_ fl ^_ e -Utf (D ?) 



which can directly be verified by differentiating the right-hand side with respect to t' . Applying 
this result to Eq. (|D.6j) we obtain 

G<(t, t) = if ^/(w - /x) J] {G R (^)r Q G A (a;) + e ih ^G R (oj + V a )r a G A (u + V^e^ 

OL 

-G R {oj + V a )r a G A (uj + V a ) 

+V a G R (ui)G R (u} + V a )r Q G A {uj)G A (u + V^e^+^^t^ 
-V a G R (oj)G R (co + V a )r a G A (uj)G A (oj + V a ) 
+V a G R {u + V^e-'^+^-^r^iu + V a )G A (u) 

-V a G R (u + V a )r a G A (co + V a )G A (uj)} . (D.8) 

Then the definition for G can be inserted into the left-hand side, and multiply accordingly with 
g-ifteff* f rom l e ft anc j w ith e lh e« l from right. Combining terms according to Eqs. (|D.lj) and (|D.2p 
we find Eq. pO]) . 
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